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Detection of turbulent thermal diffusion of particles in numerical simulations 
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The phenomenon of turbulent thermal diffusion in temperature-stratified turbulence causing a 
non-diffusive turbulent flux of inertial and non-inertial particles in the direction of the turbulent 
heat flux is found using direct numerical simulations (DNS). In simulations with and without gravity, 
this phenomenon is found to cause a peak in the particle number density around the minimum of 
the mean fluid temperature for Stokes numbers less than 1, where the Stokes number is the ratio of 
particle Stokes time to turbulent Kolmogorov time at the viscous scale. Turbulent thermal diffusion 
causes the formation of large-scale inhomogeneities in the spatial distribution of inertial particles. 
The strength of this effect is maximum for Stokes numbers around unity, and decreases again for 
larger values. The dynamics of inertial particles is studied using Lagrangian modelling in forced 
temperature-stratified turbulence, whereas non-inertial particles and the fluid are described using 
DNS in an Eulerian framework. 

PACS numbers: 47.27.tb, 47.27.T-, 47.55.Hd, 47.55.Kf 
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I. INTRODUCTION 

Transport and mixing of small particles (aerosols and 
droplets) in turbulent fluid flow is of fundamental impor- 
tance in a large variety of applications (environmental 
sciences, physics of the atmosphere and meteorology, in- 
dustrial turbulent flows and turbulent combustion; see, 
e.g., [H-O). There are also astrophysical applications, 
in particular in the context of protoplanetary accretion 
discs [1-G2. 

Numerous laboratory [1, [i"^ - [i~7j and numerical fl8l — 
\2m experiments as well as observations in atmospheric 
M H 111 HI and astrophysical [Ml, IH HI turbulent 
flows have shown different kinds of large-scale and small- 
scale long-living inhomogeneities (clusters) in the spatial 
distribution of particles. It is well known that turbu- 
lent diffusion causes destruction of large-scale inhomo- 
geneities in the spatial distributions of particles. But 
how can we explain the opposite process resulting in a 
formation of large-scale clusters of particles? 

One of the mechanisms of formation of particle inho- 
mogeneities in temperature-stratified turbulence is the 
phenomenon of turbulent thermal diffusion [29]. This 
effect consists of a turbulent non-diffusive flux of iner- 
tial particles in the direction of the turbulent heat flux, 
so that particles are accumulated in the vicinity of the 
mean temperature minimum. The particular form of 
the flow field does not play any role in this effect. It 
is a purely collective phenomenon caused by tempera- 
ture stratified turbulence resulting in a pumping effect, 
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i.e., appearance of the non-zero mean effective velocity 
of particles in the direction opposite to the mean tem- 
perature gradient. A competition between two different 
phenomena, namely the turbulent thermal diffusion and 
the turbulent diffusion determines the conditions for the 
formation of large-scale particle clusters in the vicinity 
of the mean temperature minimum. The characteristic 
scale of the particle inhomogeneity formed due to the tur- 
bulent thermal diffusion is much larger than the integral 
scale of the turbulence. Furthermore, the characteristic 
time scale of the formation of the particle inhomogene- 
ity is much longer than the characteristic turbulent time, 
i.e., this is a mean-field effect. 

The phenomenon of turbulent thermal diffusion has 
been predicted theoretically in [2^] and detected in dif- 
ferent laboratory experiments in stably and unstably 
temperature-stratified turbulence |3fj| - l33 |. This phe- 
nomenon is shown to be important for atmospheric tur- 
bulence with temperature inversions [35j and for small- 
scale particle clustering in temperature-stratified turbu- 
lence [34], but it is also expected to be significant for 
different kinds of heat exchangers, e.g., industrial boilers 
where Reynolds numbers and temperature gradients are 
large. 

In spite of the fact that turbulent thermal diffusion 
has already been found in different types of laboratory 
experiments and atmospheric flows, this effect has never 
been observed in direct numerical simulations. The main 
goal of this paper is to find turbulent thermal diffusion 
of non-inertial and inertial particles in direct numerical 
simulations (DNS). 

The paper is organized as follows. In Sect. II we dis- 
cuss the physics of the phenomenon of turbulent thermal 
diffusion. The numerical simulations for fluid, inertial 
and non-inertial particles, and the results of direct nu- 
merical simulations are described in Sect. III. Motions 
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of inertial particles are determined using a Lagrangian 
framework (Sections III-B,C,D), while non-inertial par- 
ticles are described using an Eulerian framework (Sect. 
III-E). Conclusions are drawn in Sect. IV. 



II. PHYSICS OF TURBULENT THERMAL 
DIFFUSION 

In this section we discuss the physics of the phe- 
nomenon of turbulent thermal diffusion. The number 
density n p (t, r) of particles advected in a turbulent flow 
is determined by the following equation [H, 113] 



DVr 



= 0, 



(1) 



where D is the coefficient of Brownian diffusion of parti- 
cles and U p is a random velocity of the particles which 
they acquire in a turbulent velocity field U. In order to 
study the formation of large-scale inhomogeneous parti- 
cle structures, Eq. (TTJ) for the particle number density 
is averaged over an ensemble of turbulent velocity fields. 
This yields an equation for the mean number density, 
N, of particles that has been derived using different ap- 
proaches in m, m m^m ; 

dN 



dt 
where 



V- [TV (W g + V cS ) - {D + D T ) VN] = , (2) 



V ett = -r f U p (V • Up) = -aD T V\nT 



(3) 



is the effective pumping velocity of particles due to turbu- 
lent thermal diffusion, T is the mean fluid temperature, 
a is the turbulent thermal diffusion ratio (see below), Dt 
is the turbulent diffusion coefficient, Tf is the fluid turbu- 
lent integral timescale, overbars denote ensemble averag- 
ing, Wg = T p g is the terminal fall velocity of particles, 
g is the gravitational acceleration, r p is the Stokes time 
that describes particle-fluid interactions. Equation ^ is 
written for a zero mean fluid velocity. For large Peclet 
numbers, the coefficient a is unity for non-inertial par- 
ticles, while for inertial particles a depends on r B , the 



Mach number and the fluid Reynolds number 

Already simple arguments allow us to estimate the ef- 
fective particle velocity V cS . Let us estimate the turbu- 



lent flux of particles nU p . To this end we average Eq. (fTJ) 
over an ensemble of turbulent velocity fields, and subtract 
the obtained averaged equation from Eq. (fTJ) . This yields 
an equation for the fluctuations of the particle number 
density, n = n p — N: 



^ + V- [nUp-nUp-DVn] 



-V-(NUp), (4) 



where for simplicity we consider the case of zero mean 
particle velocity, U p = 0. The left hand side of Eq. (f?]) 
has dimension n/r, where the time r can be identified 
with the fluid turbulent integral timescale Tf for large 



fluid Reynolds numbers. The fluctuations of particle 
number density can then be estimated as 



- Ti (V-Up)N - Ti (Up-V)N. 



(5) 



Multiply Eq. ([5]) by U p and average this equation over 
an ensemble of turbulent velocity fields, we arrive at a 
formula for the turbulent flux of particles: 



nUp = NV C 



D T VJV, 



(6) 



where the first term on the right hand side of Eq. ([6]) 
describes the turbulent flux of part icles due to t urbulent 
thermal diffusion: NV ctf = -Nt { U p (V ■ Up), and the 
second term on the right hand side of Eq. (J6)) determines 
the flux of particles caused by turbulent diffusion (see, 
e.g., |H): 



D T VN = rf (U p )i (U p )j VjN « r f U t U 3 VjJV, (7) 



and Dt = Tf Ui Uj is the turbulent diffusion tensor. For 
isotropic turbulence Dt = DtS^, where Dt = Tf U 2 /3. 

A. Inertial particles 

The non-diffusive mean flux of particles, NV , to- 
ward the mean temperature minimum is the main rea- 
son for the formation of large-scale inhomogeneous dis- 
tributions of inertial particles in temperature-stratified 
turbulence. Indeed, the steady-state solution of Eqs. ([2]) 
and ([3]) for the mean number density of inertial particles 
is given by 



N(z) = [f{z)]~ a exp 



[W S /D T ] dz' 



(8) 



where W g = |W g | is the modulus of the terminal fall 
velocity, T — T/Tq is the non-dimensional mean fluid 
temperature, To is the mean temperature far from the 
cooling zone, and N = N/Nq is the non-dimensional 
mean number density of particles. Here the subscripts 
represent the values at the boundary, and we neglect 
the small molecular flux of particles in Eq. © which de- 
scribes the Brownian diffusion of particles. Equation ([8|) 
implies that small particles are accumulated in the vicin- 
ity of the mean temperature minimum. 

In the following the mechanism of turbulent thermal 
diffusion for inertial particles with material density, p p , 
much larger than the fluid density, p, will be explained 
[29l | . The inertia causes particles inside the turbulent ed- 
dies to drift out to the boundary regions between eddies. 
This can be seen by considering the equation of motion 
for particles, dU p /dt = —(U p — U)/t p , where, for sim- 
plicity, the gravity force is neglected. The solution of the 
equation of motion for small Stokes time reads (4lj : 



U n = U 



T, 



dU 



0(r. 



(9) 
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For large Reynolds numbers, this yields: 

V-C/p = V-C/+^ V 2 p + 0(r p 2 ), (10) 
P 

where p is the fluid pressure. For large Peclet numbers, 
when molecular diffusion of particles in Eq. ([1]) can be 
neglected, it follows that V-C/ p oc — dlnn/di. There- 
fore, in regions with maximum fluid pressure (i.e., where 
V 2 p < 0), there is accumulation of inertial particles [i.e., 
dn/dt oc —N(T p /p)V 2 p > 0]. These regions have low 
vorticity, high strain rate, and maximum fluid pressure. 
Similarly, there is an outflow of inertial particles from 
regions with minimum fluid pressure. 

In homogeneous and isotropic turbulence without 
mean gradients of temperature, a drift from regions with 
increased concentration of particles by a turbulent flow 
is equiprobable in all directions, and the pressure (tem- 
perature) of the surrounding fluid is not correlated with 
the turbulent velocity field. There is only non-zero cor- 
relation (U ■ V)p which contributes to the energy flux, 
while Up = 0. 

In temperature-stratified turbulence, temperature and 
velocity fluctuations are correlated due to a non-zero tur- 
bulent heat flux, U 9 ^ 0, where 9 is the fluid tempera- 
ture fluctuation. Fluctuations of temperature cause pres- 
sure fluctuations, which result in fluctuations of the num- 
ber density of particles. 

Increase of the pressure of the surrounding fluid is ac- 
companied by an accumulation of particles, and the di- 
rection of the mean flux of particles coincides with that 
of the turbulent heat flux. The mean flux of particles 
is directed toward the minimum of the mean tempera- 
ture, and the particles tend to be accumulated in this 
region [2j| [H, [38[ . To demonstrate that the directions 
of the mean flux of particles and the turbulent heat flux 
coincide, we assume that the mean temperature Ti at 
point 2 is larger than the mean temperature T\ at point 
1 (see, e.g., [32||). We consider two small control volumes 
"a" and "b" located between these two points; see Ta- 
ble |U Let the direction of the local turbulent velocity in 
control volume "a" at some instant be the same as the 



TABLE I. Mechanism of particle transport n U for inertial 
particles along the direction of turbulent heat flux. 



» U9 
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b 


U px > 


U px < 


6 > 


6 < 


p > 


p < 


V- U p < 


V- Up > 


n > 


n < 


n U px > 


n U p x > 



■> nUp 



direction of the turbulent heat flux U 9 (i.e., along the x 
axis toward point 1) and let the local turbulent velocity 
in control volume "b" , at the same instant, be directed 
opposite to the turbulent heat flux (i.e., toward point 2). 

In temperature stratified turbulence, fluctuations of 
pressure p and velocity U are correlated, and regions 
with a higher level of pressure fluctuations have higher 
temperature and velocity fluctuations. The fluctuations 
of temperature 9 and pressure p in control volume "a" are 
positive because U px 9 > 0, and negative in control vol- 
ume "b" ; see Table|TJ The fluctuations of particle number 
density n are positive in control volume "a" [because par- 
ticles are locally accumulated in the vicinity of the maxi- 
mum of pressure fluctuations, dn/dt oc — N (t p /p) V 2 p > 
0], and they are negative in control volume "b" (because 
there is an outflow of particles from regions with low 
pressure) . The flux of particles n U x is positive in con- 
trol volume "a" (i.e., it is directed toward point 1), and 
it is also positive in control volume "b" (because both 
fluctuations of velocity and number density of particles 
are negative in control volume "b" ). Therefore the mean 
flux of particles n U is directed, as is the turbulent heat 
flux U9, toward point 1. This causes the formation of 
large-scale inhomogeneous structures in the spatial dis- 
tribution of inertial particles in the vicinity of the mean 
temperature minimum. 

B. Non-inertial particles 

Let us now consider non-inertial particles advected by 
the fluid flow such that the particle velocity coincides 
with the fluid velocity, U p = U. For low Mach numbers, 
using the anelastic approximation, V • U rs — U ■ Vlnp, 
the effective velocity of non-inertial particles is given by 

Vf = Ui (V • XT) = T f U~U~V 3 Inp, (11) 

where ~p is the mean fluid density. When Vp = and 
anisotropy of turbulence is weak, the effective velocity of 
non-inertial particles is: 

V cS = D T V lnp= -Z>r VlnT, (12) 

where p is the mean fluid pressure. The steady-state 
solution of Eqs. © and ([12")) for the mean number density 
of non-inertial particles is given by 

N(z) = [p(z)] DT ^ D+DT \ (13) 

where p{z) — ~p/~p~o, where po is the fluid density at the 
boundary. 

It follows from Eq. (TT2l that the non-diffusive flux of 
particles is directed toward the maximum of the fluid 
density. The physics of the accumulation of the non- 
inertial particles in the vicinity of the maximum of the 
mean fluid density (or the minimum of the mean fluid 
temperature) is quite different from that of inertial parti- 
cles and can be explained as follows (see Table HH. Let us 
assume that the mean fluid density p 2 a t point 2 is larger 
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TABLE II. Mechanism of particle transport n U for non- 
inertial particles along the direction of the mean fluid density 
gradient. 

► Vp 



a 


b 


U x >0 


U x < 


V-C7 < 


V-LT > 


n > 


n < 


nU x > 


nt/j > 



-» rat/ 



than the mean fluid density ~p 1 at point 1. Consider two 
small control volumes "a" and "b" located between these 
two points, and let the direction of the local turbulent 
velocity in control volume "a" at some instant be the 
same as the direction of the mean fluid density gradient 
V/5 (i.e., along the x axis toward point 2). Let the local 
turbulent velocity in control volume "b" be directed at 
this instant opposite to the mean fluid density gradient 
(i.e., toward point 1). 

In a fluid flow with an imposed mean temperature 
gradient (i.e., an imposed mean fluid density gradient), 
one of the sources of particle number density fluctu- 
ations, n oc — Tf N (V- U), is caused by the non-zero 
V • U fa — [7 • Vln/3 ^ [see the first term on the 
right hand side of Eq. fl5}]. Since the fluctuations of 
the fluid velocity U are positive in control volume "a" 
and negative in control volume "b", we have V- U < 
in control volume "a", and V- U > in control volume 
"b" . Therefore, the fluctuations of the particle number 
density n oc — Tf N (V- U) are positive in control volume 
"a" and negative in control volume "b" . However, the 
flux of particles n U x is positive in control volume "a" 
(i.e., it is directed to point 2), and it is also positive in 
control volume "b" (because both fluctuations of fluid 
velocity and number density of particles are negative in 
control volume "b"). Therefore, the mean flux of parti- 
cles n U is directed, as is the mean fluid density gradient 
Vp, toward point 2. This forms large-scale heteroge- 
neous structures of non-inertial particles in regions with 
a mean fluid density maximum. 



fluid density p, the fluid velocity U, and the specific en 
tropy s: 

Din p 



Dt 



= -V-U., 



^ = -- p [Vp-V.(2puS)} 



Ds 1 
T— = -V 



KVT + 2vS 2 -c P {T- T ref ), 



(14) 
(15) 
(16) 



Dt p 

where T = T exp(s/cv)(p/po) 7_1 is the fluid tempera- 
ture, 7 = cp/cy is the ratio of specific heats at constant 
pressure and constant volume, respectively, D/Dt = 
d/dt + U-'V is the advective derivative, / is the external 
forcing function, p = (cp — cy)pT is the fluid pressure, v 
is the kinematic viscosity, and K is the thermal conduc- 
tivity. The traceless rate of strain tensor is given by 



where Uij — VjiX, 



The last term in the entropy equa- 
tion (|16p determines the cooling, which causes a temper- 
ature minimum at z — 0, where 

r re£ = T + ST cxp(-z 2 /2 ( j 2 ), 

z is the vertical coordinate. 

The size of the simulation domain is L in all three 
directions. The smallest wavenumber in the domain, 
k% = 2n/L. Wavenumbers are measured in units of fci, 
and lengths in units of . The width of the cooling 
function is a = 0.5/fci. Furthermore, the strength of the 
cooling function, ST, is such that the relative contrast 
between temperature maximum and temperature mini- 
mum is ~ 1.6. In all simulations the Prandtl number 
Pr = v/K = 1. 

Turbulence in the simulation box is produced by the 
forcing function /, which is solenoidal and non- helical, 
i.e., V-/ = /- Vx/ = 0, and injects energy and mo- 
mentum perpendicular to a random wavevector whose 
direction changes every timestep, but its length is ap- 
proximately fcf (see j44()- The strength of the forcing is 
set such that the maximum Mach number is below 0.5. 
It has been checked that decreasing the Mach number 
by a factor of 5 does not have a significant effect on the 
results. 



B. Lagrangian model for inertial particles 



III. DIRECT NUMERICAL SIMULATIONS 
A. DNS model for the fluid 

In this study, equations for the fluid are solved by em- 
ploying DNS in an Eulerian framework. All numerical 
simulations have been performed using the Pencil Code 
42] (for details of the code, see Refs. [Mill)- The set 
of compressible hydrodynamic equations is solved for the 



The equation of motion of the inertial particles is 
solved numerically in a Lagrangian framework. Particles 
are treated as point particles and we consider the one- 
way coupling approximation, i.e., there is an effect of the 
fluid on the particles only, while the particles do not in- 
fluence the fluid motions. This is a good approximation 
when the spatial density of particles is much smaller than 
the fluid density. The particle equation of motion is 



dC/p 
dt 



u p -u 



(17) 
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where dX/dt = U p 



p 37T/Xrf(l-/ C ) 

is the Stokes time, m p = p p (ird 3 /6) is the particle mass 
of the spherical form, d is the particle diameter, p p is 
the particle material density, p = p v is the dynamic vis- 
cosity, / c = 0.15Re p 687 Q, and Re p = \U P - t/|d/z/ is 
the particle Reynolds number. The key parameter of the 
problem is the Stokes number St = T p /rk that is based on 
the particle Stokes time and the Kolmogorov timescale 
7"k = Tf/VRe, where n = l/u rms fcf is the turbulent inte- 
gral timescale, Re = u rms /vk^ is the fluid Reynolds num- 
ber and M rms is the root mean square fluid velocity. The 
parameter f c is used for inertial particles with Re p > 1. 
Even though f c is indeed used for calculating the particle 
forces in the simulations it is set to zero when calculating 
the Stokes number of a given particle. This is justified 
by the fact that for all particles considered here f c is al- 
ways very small and, furthermore, the Stokes number is 
only used for illustrative purposes. The ratio between 
the particle material density and the fluid mass density 
is S = Pp/po, and for all simulations S — 1000 outside 
the cooled zone. 

To demonstrate robustness of the effect of particle ac- 
cumulation in the vicinity of mean temperature minimum 
we use three different models for the dynamic viscosity 
p: (i) p — p t y / T(z) /To, where p* =const; (ii) /x=const; 
(iii) p = p(z) v, where v =const. 

For most of the simulations the number of Eulerian 
grid points are 128 3 , while the number of Lagrangian 
particles is 5 x 10 5 . It has been verified that increasing 
the number of grid points up to 512 3 or the number of 
particles up to 10 7 has no effect on the results. 

In simulations without gravity, periodic boundary con- 
ditions are used in all directions both for the fluid and for 
the particles. When gravity is taken into account in the 
simulations, particles are made elastically reflecting from 
the vertical boundaries. It has been checked that adding 
small Brownian diffusion of particles does not affect the 
results. 

In all our simulations, gravity is ignored in Eq. (|15[) for 
the fluid, but for the simulations presented in § 1111 CI it 
is included in Eq. (|17[) for particle motions. This situa- 
tion has direct applications to atmospheric turbulence 
with lower-troposphere temperature inversions, where 
the density scale height due to gravity is about 8 km, 
while the characteristic temperature inhomogeneity scale 
inside the temperature inversions is about 500-800 m, 
and the integral scale of turbulence is about 50-100 m 
(see, e.g., @). 



C. Inertial particles in the absence of gravity 

Let us first discuss the results of the numerical simula- 
tions concerning the formation of the inertial particle in- 
homogeneities in the absence of gravity. In Figs. [T] and [5] 
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FIG. 1. Upper plot: vertical profile of the mean number 
density of inertial particles for different Stokes numbers for 
simulations with kf = 5/si, Re = 240 in the case where the 
kinematic viscosity is kept constant and no gravity. Here the 
Stokes number St is based on the fluid density and temper- 
ature at the boundary. Lower plot: same as the upper plot, 
but showing the particle number density only for the smallest 
Stokes numbers together with the fluid mass density. 
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FIG. 2. The same as the upper plot of Fig. [T] but for the 
case where the dynamic viscosity is p = p,* y/T(z) /To, p* is 
kept constant, and the Stokes number St is based on the fluid 
density and temperature at the boundary. 
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we plot the vertical profile of the mean number density 
of particles for simulations with kf = 5fci, Re = 240, dif- 
ferent Stokes numbers, and different formulations for the 
dynamic viscosity. Inspection of these figures shows that, 
for a wide range of Stokes numbers, the maximum in the 
particle number density is located in the vicinity of the 
mean temperature minimum. 

In the absence of gravity (W g = 0), the parameter 
a is given by a = - In N / In T; see Eq. ©. In Fig. [3] 
the Stokes number dependence of the parameter a is 
shown for different Reynolds numbers and different forc- 
ing scales. It can be seen that the parameter a reaches 
its maximum value for St* « 0.2, while for Stokes num- 
bers St > St*, the value of a decreases monotonically. 
A critical Stokes number, St c , separating the large and 
small Stokes number regimes might be defined such that 
a > 1 for St < St c , while for St > St c the parameter a is 
decreasing steeply with increasing Stokes number. 

Since in the simulations the rms Mach number is not 
very small (around 0.2), the maximum value of the 
parameter a is only slightly larger than 1. Indeed, 
the contribution of particle inertia to the coefficient a 
for small Stokes numbers is determined by V-C7 p cx 
(r p /p) V 2 p. On the other hand, V 2 p/p cx c 2 V 2 6/T cx 



- 2 u 2 

"rms 



V z 6/T, where M 



s /c s is the Mach num- 



M 

ber, and c s is the speed of sound. The deviation of the 
parameter a from 1 is of the order of [34| 
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FIG. 3. The parameter a versus the Stokes number St 
for different Reynolds numbers and forcing wavenumbers in 
the absence of gravity. For the dashed-triple-dotted line the 
dynamic viscosity, /i, is kept constant while for the dotted 
line named "temp" the dynamic viscosity n — jj,,y / T(z) /To, 
where is kept constant. For all other lines the dynamic 
viscosity, fi — p(z)v, where the kinematic viscosity v is kept 
constant. Here the Stokes number St is based on the values 
of fluid density and temperature at the boundary. 



1 cx 



— M~ 



(18) 



In this estimate we took into account that \U"V 2 9\ ~ 
i~ 2 \U9\ and the turbulent heat flux is TJO = -k t VT, 
where n T ~ Tf u 2 ms and £ = Tf u rnls . 

Therefore, decreasing the Mach number causes a to 
increase due to particle inertia. For instance, in atmo- 
spheric turbulence the Mach number is of the order of 
(0.1-3) xl0~ 3 and a is about 10 35]. On the other hand, 
the observed effect of accumulation of inertial particles 
occurs only for turbulent flows with large fluid Reynolds 
number. Decreasing the Mach number in the Pencil 
Code will either result in a shorter time step or a lower 
Reynolds number. This implies that we cannot easily 
reach at the same time small Mach numbers and large 
fluid Reynolds number. Therefore, the observations of 
the accumulation of inertial particles in the parameter 
range a 3> 1 is not currently possible with the Pencil 
Code. 

The parameter a for the heavy particles decreases 
strongly with r p cx St (see Fig. [3] in the large Stokes 
number regime). Furthermore, it is expected that, in 
the large Stokes number regime, a is independent of Re 
for a particle as long as u rms and kf are kept constant. 
This is due to the fact that particles with t p ~ Tf are 
affected primarily by the largest turbulent eddies in the 
flow. This is because the smaller turbulent eddies have 
turnover times much shorter than the particle's Stokes 
time, and they cannot accelerate the particles. Since the 



Stokes number is based on the Kolmogorov time t^, this 
implies that St c cx %/Re, which is indeed what is found 
in Fig. H 

When experiencing a high fluid density, a particle in 
the large Stokes number regime will more easily be ac- 
celerated by turbulence if the kinematic viscosity v is 
constant than if the dynamic viscosity /i = pv is con- 
stant. This is due to the fact that for constant v a high 
fluid density is associated with a smaller Stokes number, 
while for constant \x the Stokes number is independent of 
the fluid density. Due to this it is expected that parti- 
cles with large Stokes numbers tend to be more depleted 
from the high density regions when v is constant than 
when fj, is constant. This explains why the simulation 
shown in Fig. [3] with constant fi has a much shallower 
fall-off with increasing Stokes number in the large Stokes 
number regime than the simulations with constant v. 

If the integral scale of turbulence is close to the scale 
of the thermal cooling layer, a particle trapped inside 
a turbulent eddy might travel across the whole cooling 
zone during one eddy turnover time. This will effectively 
smear out the peak of the particle number density, and 
consequently also decrease the parameter a. For simula- 
tions with kf — 5fci, the integral scale is comparable to 
the cooling scale, and the peak of a for kf = 5k\ is lower 
than for simulations with kf = 15fci (see Fig. This 
trend is expected to continue for yet larger kf. 

In addition, Fig. [3] shows that the difference in the pa- 
rameter a determined for different formulations of the 
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TABLE III. The coefficient a, the relative maximum of the 
particle number density and profile width for simulations with 
kf — 5fei, Re = 240 and for different formulations of the 
dynamic viscosity: (i) in the model /i the dynamic viscosity 
is fi = /x* \YT(z)/To, where p* is constant; (ii) in the model 
/2 the dynamic viscosity is constant; (iii) in the model fs the 
dynamic viscosity p — p(z) v, where the kinematic viscosity v 
is constant. The profile width Lm is given at the level when 
the number density is decreased in e times. 







a 




N max /N Q 


Profile width 


St 


h 


h 


h 


h 


h h 


h 


h 


h 


0.1 


1.00 


1.01 


1.00 


1.55 


1.66 1.63 


0.88 


0.80 


0.79 


1 


0.91 


0.98 


0.95 


1.49 


1.62 1.60 


0.86 


0.80 


0.79 


10 


0.80 


0.75 


0.59 


1.42 


1.44 1.35 


0.90 


0.82 


0.81 


14 


0.76 


0.68 


0.45 


1.39 


1.39 1.25 


0.90 


0.83 


0.82 



dynamic viscosity p, is very small for St -C 1, and it is 
small for St ~ 1. The real difference in the parameter 
a is only observed for large Stokes numbers, St > 10 
(see TABLE In this case, the largest value of a oc- 
curs for the dynamic viscosity p = p*^/T(z)/Tq, where 
p* =const. On the other hand, the difference in relative 
maximum of the particle number density and the profile 
width for different formulations of dynamic viscosity is 
very minor (see TABLE IIII[) . Also scale-separation in- 
creases effect of the particle accumulation for St > 1 [see 
Fig. [3]. However, using different parameters and different 
formulations of the dynamic viscosity, we always observe 
the effect of particle accumulation in the vicinity of the 
mean temperature minimum. This is an indication of the 
robustness of this phenomenon. 



D. Inertial particles with gravity 

For heavy inertial particles with large Stokes numbers, 
gravity plays a crucial role, and must be included in the 
simulations. In Fig. 01 the mean particle number density 
is shown as a function of vertical position z for simula- 
tions with St = 1, k{ = 5ki, and different gravitational 
accelerations. 

In the following, gravity is measured in the units of 
go = Dt/^t^L), so that g/ga — L/L g for a particle 
of Stokes number unity. Here, L g — D^/Wg is the 
characteristic scale of the mean particle number den- 
sity variations due to the gravity for an isothermal case, 
£>t = u rms /3fcf is the turbulent diffusion coefficient, and 
L is the height of the box. 

As the particle sedimentation velocity is increased, the 
particle number density profile is more and more tilted, 
as expected. For large particle sedimentation velocity 
(9/90 — 200/3), almost all particles have accumulated at 
the lower wall of the box (see Fig. @|. 

In Fig. [SJ the vertical profile of the particle number 
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FIG. 4. Vertical profile of the mean number density for dif- 
ferent gravitational accelerations, St = 1 and kt = 5fci, where 
the kinematic viscosity v is kept constant. 



density is shown for .9/50 = 2/3 and different Stokes num- 
bers. By fitting these results with Eq. ©, the parameter 
a is found to be about 1 for the three smallest Stokes 
numbers in Fig. [5J while for St = 3.5 and St = 14 a 
best fit of a is found around 0.8 and 0.5, respectively. 
For large Stokes numbers the maximum particle number 
density is found near the bottom wall of the box due to 
the large sedimentation velocity. 



E. Non-inertial particles 

For comparison with the numerical simulations per- 
formed for inertial particles in a Lagrangian framework, 
we describe in this section numerical simulations for non- 
inertial particles, where the equations for both, fluid and 
particles are solved by employing DNS in an Eulerian 
framework. In particular, we now solve Eq. ([T]) for the 
number density of non-inertial particles n p with U p = U , 
and the fluid density p, the fluid velocity U, and the spe- 
cific entropy s, using again the Pencil Code For 
these simulations the resolution is 128 3 . For the fluid we 
apply the same conditions as was described in previous 
subsections. In particular, the periodic boundary condi- 
tions are used in these simulations in three directions for 
Eqs. © and 

Non-inertial particles are characterized by the follow- 
ing dimensionless parameters: Pe = u rms /Dk{ is the 
Peclet number and Sc = v/D is the Schmidt number. 
The results of the simulations, shown in Fig. [51 demon- 
strate the accumulation of non-inertial particles in the 
vicinity of the maximum of the mean fluid density (or 
the minimum of the mean fluid temperature). This is in 
agreement with the theoretical predictions |29j . 

The results shown in Fig. [^demonstrate that the pro- 
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St= 0.05 



an Eulerian framework. The phenomenon of turbulent 
thermal diffusion has been studied for different Stokes 
and fluid Reynolds numbers, Peclet numbers as well as 
different forcing scales of the turbulence. Furthermore, 
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FIG. 5. Vertical profile of the mean number density (solid 
line) for g/go = 2/3, k{ = 5fci and different Stokes numbers 
which are based on the fluid density and temperature at the 
boundary. Here the kinematic viscosity v is kept constant. 
The dotted line represent the isothermal reference case. 



files of the mean fluid density and the mean particle 
number density are similar in numerical simulations with 
large Peclet and Reynolds numbers. However, when the 
Peclet number, Pe = Re Sc is not large, the molecular 
diffusion D becomes important, and the profiles of the 
mean fluid density and the mean particle number den- 
sity are, according to Eq. (|I3[) . no longer similar. 




FIG. 6. Vertical profiles of the mean number density (upper 
panel) and the mean fluid density (lower panel) of non-inertial 
particles for simulations with k{ = 5fci and different Reynolds 
and Schmidt numbers. 



the effect of gravity has been included in simulations. 

In all simulations, with different parameters and dif- 
ferent formulations of the dynamic viscosity, we always 
observe the effect of particle accumulation in the vicin- 
ity of the mean temperature minimum due to turbulent 
thermal diffusion for St < 1. This effect is robust and 
the results of the numerical simulations are in agreement 
with theoretical studies [H, [38| , laboratory experiments 
[3ll [33l [13] and atmospheric observations [35|. When 
St > I , this effect is decreasing with Stokes number. 



IV. CONCLUSIONS 
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